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Abstract Phylogeographic studies of Eremias lizards (Lacertidae) in East Asia have been limited, and the impact of 
major climatic events on their population dynamics remains poorly known. This study aimed to investigate population 
histories and refugia during the Last Glacial Maximum of two sympatric Eremias lizards (E. argus and E. brenchleyi) 
inhabiting northern China. We sequenced partial mitochondrial DNA from the ND4 gene for 128 individuals of Е. 
argus from nine localities, and 46 individuals of E. brenchleyi from five localities. Forty-four ND4 haplotypes were 
determined from Е. argus samples, and 33 from E. brenchleyi samples. Population expansion events began about 0.0044 
Ma in Е. argus, and 0.031 Ma in E. brenchleyi. The demographic history of E. brenchleyi indicates a long-lasting 
population decline since the most recent common ancestor, while that of E. argus indicates a continuous population 
growth. Among-population structure was significant in both species, and there were multiple refugia across their range. 
Intermittent gene flow occurred among expanded populations across multiple refugia during warmer phases of the 
glacial period, and this may explain why the effective population size has remained relatively stable in E. brenchleyi and 
grown in E. argus. 
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1. Introduction refugia within its range, the effective population size of 
the species might remain stable or grow during the glacial 
Historic climatic change (e.g., climatic oscillations during period due to intermittent gene flow between refugia 
during warmer periods (interstadial) (Li et al., 2009; 


Ruiz-González et al., 2013). Climatic changes during 


the Pleistocene Epoch) affected the historical demography 
of extant species (Hewitt, 2000, 2004). Glacial and 


interglacial periods consisted of a series of alternating Pleistocene cycles have played a key role in shaping 


genetic diversity of species and left genetic signatures 
in present populations (Hewitt, 1996, 1999, 2000, 2004; 
Avise, 2000) that can help identify glacial refugia, 


warming periods of limited duration (interstadials) and 
short-term cooling periods (stadials) (Martrat et al., 


РЛ ЕЕ ЕЕЕ QUOD онаа post-glacial migration routes and population dynamics 


decrease during the glacial period and increase during 
the interglacial period. However, if there were multiple 
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(Newton et al., 1999; Petit et al., 2005; Li et al., 2009; 
Tian et al., 2009; Ruiz-González et al., 2013). Relative 
to birds and mammals, reptiles are less able to disperse 
and are therefore more sensitive to climatic fluctuations, 
making them more suitable biogeographical indicators 
(Camargo et al., 2010; Malhotra et al., 2011). 
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Lacertid lizards of the genus Eremias are widely 
distributed in Southeast Europe, West and Central Asia, 
and Mongolia, northern China, Korea and Russia (regions 
of Lake Baikal) in East Asia. Of some 30 extant species, 
eight can be found in China, six of which (Е. argus, 
E. arguta, E. brenchleyi, E. grammica, E. velox and E. 
vermiculata) are oviparous, and two (E. multiocellata and 
E. przewalskii) are viviparous (Zhao, 1999). The ancestor 
of Eremias lizards has been hypothesized to originate 
from Europe and dispersed northwards to Central Asia 
through Africa during the late Miocene 8-10 million 
years ago (Ma) (Ға, 1998; Arnold et al., 2007; Mayer and 
Pavlicev, 2007; Hipsley et al., 2009). According to this 
hypothesis, Eremias lizards found in East Asia should 
spread from west (Central Asia) to east, not following 
westward, southward or northward adaptive radiation of 
species found in other parts of the world (Forcina et al., 
2012; Mori et al., 2013; Zink et al., 2013). During the Last 
Glacial Maximum (LGM), multiple glacial refugia were 
available to East Asian species throughout their ranges 
(Qian and Ricklefs, 2000; Li et al., 2009; Tian et al., 
2009; Bai et al., 2010) primarily due to a relatively mild 
Pleistocene climate in that region (Weaver et al., 1998; 
Ju et al., 2007). However, as phylogeographic studies of 
Eremias lizards in East Asia have been limited, the impact 
of major climatic events on their population dynamics 
remains poorly known. 

Here, we study two oviparous Eremias lizards (E. 
argus and E. brenchleyi) sympatric throughout the 
distribution of E. brenchleyi in the central and northern 
parts of China (Chen, 1991). The two species can coexist 
primarily because they use different microhabitats (Е. 
argus mainly inhabits plains and hills, while E. brenchleyi 
inhabits hills only) not because they differ greatly 
along the resource axes of diet and time (Zhao, 1999). 
Both species are relatively common and have more 
intact population genetic structure than other species or 
populations that have become threatened or endangered in 
China (Elderkin et al., 2007; Lin et al., 2010; Zhao et al., 
2011), thus providing ideal model systems in which to 
investigate population histories, glacial refugia and post- 
glacial migration routes. Here, we used partial mtDNA 
sequences from the gene ND4 to address three questions 
regarding population dynamics following LGM in these 
two species. (1) Did they survive in multiple refugia 
across their species range or in a single refuge during the 
LGM? (2) What was the pattern of recolonization routes? 
(3) How did effective population size vary during the 
LGM? 


2. Materials and Methods 


2.1 Sampling We sequenced the mitochondrial ND4 gene 
from 128 individuals of E. argus from nine sites (Figure 
1 A, Table 1) and 46 individuals of Е. brenchleyi from 
five sites (Figure 1 B, Table 1). The most distal 2-3 mm 
of the tail tip from each lizard was stored in 95% ethanol. 
Lizards were released at their site of capture immediately 
after tissue sampling. The influence of tail breaks at the 
distal portion of the tail often is not significant in lizard 
species with caudal autotomy, including Eremias lizards 
(Lin and Ji, 2005; Lin et al., 2006; Zhao et al., 2008; Sun 
et al., 2009; Ding et al., 2012). Voucher specimens were 
deposited at Nanjing Normal University and assigned 
voucher numbers identified by locality-haplotype 
numbers. The sampling code, geographic information and 
genetic analysis are given in Appendix A, Appendix B 
and Table 1. We used all sequences from this study and 
three sequences from two lacertid lizards (Lacerta viridis 
viridis: МС 008328; Timon lepidus: 00902321 and 
DQ902322) as ingroups, and two sequences (DQ150379 
and DQ150376) from Psammodromus algirus (Hipsley 
et al., 2009) as outgroups for the phylogenetic analysis 
based on ND4 haplotypes. 


2.2 Molecular data Genomic DNA was isolated from 
tail samples using standard phenol-chloroform extraction 
methods (Palumbi, 1996). Fragments of complete ND4 
gene from the mitochondrial genome were amplified 
with the primers ND4 (5'-CAC CTA TGA CTA CCA 
AAA ССТ САТ GTA GAA GC-3') and leu (5'-CAT 
TAC TTT TAC TTG GAT TTG CAC CA-3') (Busack 
et al., 2005). PCRs were performed in 100 ul volumes, 
using a hot start method; PCR cycling parameters were 
94 °C for 7 min; 40 cycles of 94 °C for 40 s, 46 °C for 
30 s, 72 °C for 1 min; and a final elongation step of 72 
°C for 7 min. PCR products were then purified with 
a spin column containing Sepacry S-400 (Amersham 
Bioscience AB, Uppsala, Sweden) and sequenced using 
both forward and reverse primers on an ABI-PRISM'™ 
310 Genetic Analyzer (Applied Biosystems Information, 
USA). Sequences were edited and aligned manually 
using BIOEDIT version 7.0.9.0 (Hall, 1999), and were 
translated to amino acids with Mega 6.0 to verify if a 
functional mitochondrial DNA sequence was obtained 
and that nuclear pseudogenes were not amplified (Tamura 
et al., 2013). All sequences were deposited in GenBank 
under accession numbers JN561171—JN561247. 


2.3 Phylogenetic analysis Phyolgenetic analysis were 
carried out by Bayesian inference (BI) using MrBayes 
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Figure 1 Map of northern China showing localities where individual E. argus (A) and E. brenchleyi (B) were collected. See text and Table 
1 for sample locality abbreviations. Numbers in parentheses are sample sizes. Small circles represent sampling localities, and pie charts 
represent haplotype group frequencies within each locality. The blue, dotted lines outline the species range in China. 


3.1.2 (Ronquist and Huelsenbeck, 2003) and MrModeltest 
2.3 (Nylander, 2004). The best-fit model (GTR + I) for 
15 codon position, (HKY + G) for 2nd codon position, 
and (GTR + G) for 3rd codon position were selected by 
Akaike Information Criterion (AIC) in MrModeltest 2.3. 
Four Markov Chains Monte Carlo (MCMC) samples 
were run for 5 x 10° iterations. Two independent runs 
were carried out to allow additional confirmation of 
the convergence of MCMC runs. Trees were sampled 


every 500 iterations, providing 2 x 10* samples from the 
two runs. The first 4% of the iterations were discarded 
as burn-in, leaving 1.96 x 10* samples to estimate the 
consensus tree and the Bayesian posterior probabilities. 
Haplotypes with lineages were evaluated using median- 
joining networks (MJNs) in NETWORK 4.5.1.0 (Bandelt 
et al., 1999). 


2.4 Divergence times We used Bayesian MCMC method 
to evaluated divergence times and evolutionary rates. For 
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Table 1 Descriptive statistics by sampling locality for the two lizards studied. NH: number of haplotypes; LSH: locality-specific haplotypes; 
HD: haplotype diversity; ND: nucleotide diversity; D: Tajima's D; Fu’s Ғо statistics of Fu’s Е; test (1997); OMMD: observed modality 
of mismatch distribution; SSD: the sum of square deviation between the observed and simulated mismatch distributions; Raggedness: 
raggedness index. 


Samplin 
Species locality Latitude Longitude п NH LSH HD (/) ND (л) D Fu's Е; OMMD SSD Raggedness 
(abbreviation) 
E. argus 128 44 0.952 0.025 -0.59 —2.83 Multimodal 0.00725 0.01277* 
Harbin (HRB) 45.7 126.6 22 7 4 0.541 0.01 -097 Sell Multimodal 0.34002*** 0.12777 
Xinghe (XH) 41.1 113.9 12 4 2 0.697 0.004 -1.56 1.48 Multimodal 0.06753 0.17034 
Etuoke (ETK) 39.5 107.5 11 8 8 0.891 0.01 -2.08*** —0.59 Multimodal 0.01838 0.04298 
Yangquan (YQ) 37.9 113.6 12 6 6 0.818 0.009 1.43 1.4 Multimodal 0.13748* 0.21212 
Handan (HD) 36.6 114.5 11 6 6 0.873 0.0142 0.52 2.5 Multimodal 0.11861** 0.13058 
Gonghe (GH) 36.4 100.5 12 2 0 0.46 0.0008 1.49 12 Unimodal 0.03392 0.30579 
Jiyuan (JY) 35.1 112.6 12 7 7 0.909 0.008 -1.64* 0.1 Multimodal 0.28364** 0.04454 
Chang'an (CA) 34 108.9 23 6 4 0.573 0.017 -0.44 8.38 Multimodal 0.38288*** 0.20803 
Chuzhou (CZ) 323 118.3 13 3 3 0.603 0.018 2.87 11.06 Unimodal 0.27213 0.45348 
E. brenchleyi 46 33 0.984 0.059 1.69 —0.88 Multimodal 0.02235 0.00581 
Yangquan (YQ) 37.9 113.6 11 8 8 0.946 0.008 0.04 -1.38 Multimodal 0.01076 0.0281 
Handan (HD) 36.6 114.5 9 6 6 0.889 0.018 1.6 1.96 Multimodal 0.09925** 0.18827 
Taishan (TS) 36.3 117.1 3. 3. 3 1 0.004 0 -0.34 Unimodal 0.23127 0.66667 
Jiyuan (JY) 35.1 112.6 11 9 9 0.946 0.004 0.78 —5.52***  Unimodal 0.00411 0.0519 
Xuzhou (XZ) 34.3 117.2 12 7 7 0.84 20.003 -0.66 —2.68* Unimodal 0.06078 0.24174 


*: P«0.05; **: Р < 0.01; ***: P < 0.001. 


E. argus, we set the mean of the prior distribution for the 
time separating the ingroup root from the present (rttm) 
to 45.4 Ma estimated based on the age (44.1—46.7 Ma 
with a mean of 45.4 Ma; Appendix A, node 5 (Hipsley 
et al., 2009) of the split between Eremias and Timon. The 
standard deviation of the prior distribution for the time 
separating the ingroup root from the present (rttmsd) was 
set at one-half of rttm. The mean and standard deviation 
of prior distribution for rate at root node (rtrate and 
rtratesd) were set to 0.008 2 substitutions per site per 
1 Ma, with the rtrate value calculated by dividing the 
median distance between the ingroup root and the tips 
by rttm; the mean and standard deviation of prior for 
brownian motion constant nu (brownmean and brownsd) 
both were set to 0.022, so that rttm х brownmean = 1, 
following recommendations accompanying the software. 
АП remaining parameters were set to default values. 
In this study, the fossil calibration point was defined to 
the most recent common ancestor of Е. argus according 
to the Middle Pleistocene (0.438—0.548 Ma) found in 
Chinese fossil layers (Li et al., 2004) (node N3 in Figure 
2). The Bayesian analysis was run for 5 x 10^ generations, 
with a sample frequency of 100 after a burn-in of 5 x 10* 
generations. 

To estimate intraspecies divergence time of Е. 
brenchleyi, we set rttm to posterior distribution of the 
divergence time of Eremias and Timon, which came from 
the above results of E. argus intraspecies divergence-time 


estimation rttmsd was set to equal rttm. The rtrate value 
was calculated by dividing the median distance between 
the ingroup root and the tips by rttm and the rtratesd 
value was set to equal rtrate. We set all the remaining 
parameters to default values. We used two calibration 
points in these estimates (N15 and N16 in Figure 3). 
The N16 calibration point was based on the fossil record 
of ancestor of E. argus from the Middle Pleistocene 
(0.438—0.548 Ma) (Li et al., 2004). The lower bound 
for the age of node N16 was 0.438 Ma, and the upper 
bound for the age of node N16 was 0.548 Ma. N15 was 
not a fossil calibration but a secondary calibration point 
from the above results of the 9596 confidence interval of 
E. argus intraspecies divergence-time estimation. The 
Bayesian analysis was run for 5 x 10* generations, with 
a sample frequency of 100 after a burn-in of 5 x 104 
generations. 


2.5 Population genetic analyses We followed Zhao 
et al. (2011) to calculate haplotype diversity (Л), mean 
nucleotide diversity (z), mismatch distributions, Fu's 
Fs and Tajima's D for each sampling site and for each 
species as a whole (1000 replicates). Nucleotide diversity 
and haplotype diversity were tested for a significant 
correlation with latitude using linear regression. We 
performed mismatch-distribution tests for goodness-of- 
fit against a null model of sudden population expansion 
calculated using SSD between observed and expected 
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Figure 2 The partitioned Bayesian phylogenetic tree based оп ND4 haplotypes for Е. argus. Nodes NO to N12 denote internal nodes of 
interest. Numbers above the branches are Bayesian posterior probabilities. Taxa are haplotypes; all haplotype designations are listed in 
Appendix C, followed by the sampling localities and numbers of individuals from each locality having that haplotype [e.g., JY(1)]. 


data, and Harpending's raggedness Index (Fu, 1997; 
Excoffier et al., 2005). To further estimate past population 
dynamics through time, we constructed BSP implemented 
in BEAST 1.6.1 (Drummond and Rambaut, 2007) with a 
relaxed molecular clock method to depict the change of 
the two species' effective population size since the time 
to the most recent common ancestor (TMRCA) of the 
sampled ND4 haplotypes respectively. For each BSP, the 
model of molecular evolution was assessed for the ND4 
sequence by AIC in MrModeltest 2.3 (Nylander, 2004). 


We selected the best-fit models (GTR + G) and (GTR + 
I + G) for E. argus and Е. brenchleyi, respectively. The 
only one calibration point based on the fossil record 
from the Middle Pleistocene (0.438—0.548 Ma) (Li et al., 
2004) was assigned to the most recent common ancestor 
(MRCA) of E. argus and for this species we assumed 
a normal prior distribution centered at 0.493 Ma with a 
standard deviation of 0.028 Ma for the MRCA. For E. 
brenchleyi, the evolutionary rate of ND4 was set to the 
ucld.mean of BSP of Е. argus above because no fossil 
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Figure 3 The partitioned Bayesian phylogenetic tree based on ND4 haplotypes for E. brenchleyi. Nodes from N13 to N27 denote internal 
nodes of interest. Numbers above the branches are Bayesian posterior probabilities. Taxa are haplotypes; all haplotype designations are listed 
in Appendix C, followed by the sampling localities and numbers of individuals from each locality having that haplotype [e.g., JY(1)]. 


calibration points are available. MCMC sampling was 
run for 3 x 10” generations, with a sample frequency of 
100 after a burn-in of 3 x 10° generations. The effective 
sample size of each parameter was measured and 
demographic plots for the two species were visualized 
in Tracer version 1.5 (Rambaut and Drummond, 2007). 
We also made the diagram of oxygen isotopic (8 ^O) 
record during the past 0.78 Ma as proxy for global ice 
distribution for each species (Imbrie et al., 1984). 

We computed the pairwise ®,, value for each 
combination of populations using the Kimura two- 


parameter model of nucleotide substitution (Kimura, 
1981), which takes into account haplotype frequencies 
and the genetic distance between haplotypes in 
ARLEQUIN. The significance of ®,, values was 
evaluated using 1000 permutations of the data and set 
at P = 0.05. We performed a Mantel test to compare 
the pairwise population structure estimates to a matrix 
of geographical distance. We used ARLEQUIN also to 
perform a two-level hierarchical analysis of molecular 
variance (AMOVA) (Excoffier et al., 1992) to test for 
structure among populations: within populations, and 
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among populations (Ф). 

One additional procedure was used to characterize 
major patterns of genetic structure across the ranges 
sampled for both species. We applied the SAMOVA 
procedure (Dupanloup et al., 2002) to identify groups of 
collection localities that were maximally differentiated 
based on ND4, and to infer genetic barriers between 
these groups. We used SAMOVA 1.0 to perform analyses 
based on 100 simulated annealing steps and examined 
the number of groups (k) that maximized genetic 
differentiation among groups (Фет). The analysis was 
run fork =2 to k 7 n (n 
brenchleyi) and the significance level of fixation indices 


8 for E. argus; n — 4 for E. 


was evaluated by 1000 permutations. 
3. Results 


3.1 Matrilineal patterns We obtained 44 Е. argus 
haplotypes (Appendix А) and 33 Е. brenchleyi haplotypes 
(Appendix B). The ratio of haplotypes relative to the total 
number of individuals sampled was 0.34 for E. argus and 
0.71 for E. brenchleyi. Neither mean nucleotide diversity 
(Е = 0.69, P = 0.42) nor the number of haplotypes (7; 
— 0.96, P — 0.35) per sampling location differed between 
the two species (Table 1). The number of haplotypes, the 
number of unique haplotypes, haplotype diversity and 
nucleotide diversity for single sampling localities were 
not correlated with latitude in either species (all P > 0.26). 
All E. argus haplotypes formed a single well-supported 
clade, which could be subdivided into three clades (EA I, 
EA II and EA III) although the phylogenetic relationships 
among these clades were not fully resolved (PP « 0.50) 
(Figure 2). Clade EA I included a single haplotype 
(EA44) from the Chuzhou (CZ) population. Clade EA II 
was poorly supported (PP — 0.66) and further subdivided 
into five clades: one poorly supported clade (EA IIa, PP — 
0.63) and four well-supported clades (EA IIb, PP — 1.00; 
EA Пс, PP = 0.98; EA Па, PP = 1.00; EA Пе, PP = 1.00). 
These five clades were mostly admixed in many local 
populations and had no geographical specificity. Clade 
EA IIa was further subdivided into two poorly supported 
clades (EA Hal, PP = 0.71; EA Па2, PP = 0.88) that were 
mostly admixed in the Handan (HD) population and had 
no geographical specificity. Clade EA III included two 
haplotypes from the Harbin (HRB) population. 

АП E. brenchleyi haplotypes formed a single well- 
supported clade, which could be subdivided into two 
well-supported clades [north (NYR) and south (SYR) of 
the Yellow River, both PP — 1.00] (Figure 3). The NYR 
clade included all haplotypes from the Yangquan (YQ), 
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HD and Jiyuan (JY) populations in NYR, and the SYR 
clade included all haplotypes from the Xuzhou (XZ) and 
Taishan (TS) populations in SYR. The division of the 
two clades coincided with the Yellow River (Figure 1). 
There was a clear population-genetic structure within 
both clades. The NYR clade could be subdivided into 
three major clades (HD4, PP = 1.00; JY, PP = 0.99; 
YQ, PP = 0.81) and three minor clades (HD1, HD2 and 
HD3) of which each included a single haplotype from 
the HD population east of the Taihang Mts (36?—40?N, 
112°-115°E). The HD4 clade included haplotypes 
all from the HD population; the JY clade included all 
haplotypes from the JY population also in east of the 
Taihang Mts; the YQ clade included all haplotypes 
from the YQ population west of the Taihang Mts. The 
HD1, HD2, HD3, HD4 and JY clades formed one group 
including all haplotypes from two local populations (HD 
and JY) east of the Taihang Mts; the YQ clade formed 
the other group including all haplotypes from the YQ 
population west of the Taihang Mts. The geographical 
division of these two groups coincided with the Taihang 
Mts (Figure 1). The SYR clade could be subdivided 
into two clades with strong nodal support (TS clade and 
XZ clade, both PP = 1.00). The TS clade included all 
haplotypes from the TS population, while the XZ clade 
included all haplotypes from the XZ population. 

For Е. argus, results of the Median-joining network 
(MJN) suggests little or no association between 
haplotypes and geography (Figure 4 A, Appendix A). 
We detected related haplotypes that formed eight groups 
(EA I, EA IIal, EA Па2, EA IIb, EA Пс, EA IId, ЕА Пе 
and EA III) in the Bayesian phylogenetic tree (Figure 
2, Figure 4 A). Clade EA IIa2 and EA IIc co-occurred 
in Xinghe (XH) and YQ. Clade EA Hal and EA IIb co- 
occurred іп HD, JY and Chang’an (СА). Clade EA Па1 
and EA Па co-occurred in HRB and СА. Сізде EA Па1, 
EA IIb, EA Па and EA Пе co-occurred іп CA. Clade EA 
Hal, EA Пс, EA Па and EA III co-occurred in HRB. For 
E. brenchleyi, however, we detected networks of related 
haplotypes that were clustered into two distinct areas 
(SYR and NYR) respectively corresponding to the SYR 
clade and NYR clade in the Bayesian phylogenetic tree 
(Figure 3, Figure 4 B, Appendix B). We further detected 
networks of related SYR haplotypes that were clustered 
into two distinct areas (TS and XZ), and detected 
networks of related NYR haplotypes that were clustered 
into two distinct areas [YQ and (HD1, HD2, HD3, HD4, 
and JY)]. These four areas [((HD1, HD2, HD3, HD4, and 
JY), YQ, TS, and XZ] correspond to clades HD1, HD2, 
HD3, HD4, and JY, clade YQ, clade TS, and clade XZ, 
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Figure 4 Median-joining networks of mtDNA haplotypes at the ND4 gene in E. argus (A) and E. brenchleyi (B). Numbers adjacent to open 
circles (network nodes) refer to haplotypes listed in Appendices A and B. Haplotype group color for E. argus is consistent with Figure 1 A, 
and haplotype group color for Е. brenchleyi is consistent with Figure 1 B. 


respectively (Figure 3). 


3.2 Divergence time Table 2 gives posterior distributions 
of divergence time at each node and 95% credibility 
intervals based on phylogenies. For Е. argus, the most 
recent common ancestor (TMRCA) of clade EA Hal 
(node N7), EA Па2 (node №8), EA IIb (node N9), EA Пс 
(node N10), EA IId (node N11), EA IIe (node N12), and 
EA III (node N5) were estimated at 0.158 + 0.063, 0.178 
+ 0.064, 0.205 + 0.087, 0.220 + 0.082, 0.144 + 0.082, 


0.184 + 0.082 and 0.100 + 0.067 Ma, respectively. For Е. 
brenchleyi, the divergence time of the NYR clade from 
the SYR clade (node N15) was dated to 3.208 + 0.743 
Ma. Within the NYR clade, the divergence time of clades 
HD1, HD2, HD3, HD4 and JY from clade YQ (node N23) 
was dated to 0.939 + 0.418 Ma. Within the SYR clade, 
the divergence time of the TS clade from the XZ clade 
(node N18) was dated to 2.095 + 0.703 Ma. The TMRCA 
of clades JY and HD (node N17) was estimated at 1.798 
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Table 2 The posterior distribution of divergence times with 95% 
confidence intervals estimated by Bayesian molecular dating. DT: 
divergence times (Ma); SD: standard deviation; CI: 9596 credibility 
interval. 


Species Node DT SD for DT CI for DT 

E. argus NO 9.195 3.13 4.420-16.562 
NI 4.809 1.427 2.605-8.141 
N2 2.292 0.878 0.976—4.367 
N3 0.498 0.031 0.441—0.546 
N4 0.376 0.076 0.225-0.512 
N5 0. 0.067 0.006—0.256 
N6 0.255 0.077 0.122-0.420 
N7 0.158 0.063 0.056—0.303 
N8 0.178 0.064 0.073-0.322 
N9 0.205 0.087 0.056—0.391 
NIO 0.22 0.082 0.074—0.395 
N11 0.144 0.082 0.024—0.336 
N12 0.184 0.082 0.055-0.372 

E. brenchleyi N13 5.057 1.482 2.646—8.532 
N14 3.86 0.87 2.163—5.561 
NI5 3.208 0.743 1.648-4.309 
N16 0.495 0.032 0.441—0.545 
N17 1.798 0.615 0.700-3.035 
N18 2.095 0.703 0.817—3.482 
NI9 1.504 0.55 0.539-2.647 
N20 1.29 0.504 0.442-2.378 
N21 1.1 0.461 0.354-2.117 
N22 0.735 0.394 0.167—1.666 
N23 0.939 0.418 0.283-1.884 
N24 0.757 0.369 0.207-1.622 
N25 0.575 0.355 0.089—1.429 
N26 0.404 0.324 0.044—1.251 
N27 0.262 0.264 0.012-0.996 


+ 0.615 Ma; the TMRCA of clades YQ (node N24), XZ 
(node N27) and TS (node N26) were estimated at 0.757 + 
0.369, 0.262 + 0.264 and 0.404 + 0.324 Ma, respectively. 


3.3 Population-genetic analysis For Е. argus, mismatch 
distributions of the Gonghe (GH) and CZ populations 
exhibit relatively unimodal patterns with non-significant 
raggedness values, matching the expected distribution 
of an expanding population well, and supporting that 
these two populations have undergone rapid demographic 
expansion. Mismatch distributions of the remaining 
seven populations show multimodal distributions. 
Significant P-values for the sum of squares deviations 
(SSD), high raggedness indexes, positive Fu’s F, and 
Tajima's D values indicate no expansion in the YQ and 
HD populations, and negative Tajima's D values suggest 
recent expansion in the HRB, XH, Etuoke (ETK), JY and 
CA populations. 

For Е. brenchleyi, mismatch distributions of the TS, JY 
and XZ populations exhibit relatively unimodal patterns 
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with on-significant raggedness values, fitting the expected 
distribution of an expanding population (non-significant 
SSD) and supporting that these three populations have 
undergone rapid demographic expansion. Mismatch 
distributions of the YO and HD populations show 
multimodal distributions. The significant P-value for 
SSD, high raggedness index, positive Fu’s F; and 
Tajima's D values indicate no population expansion in 
the HD population, and negative Fu’s F; value, non- 
significant SSD and low raggedness index indicate recent 
population expansion in the YQ population. 

For all populations of E. argus, the multimodal 
mismatch distribution, non-significant P-value for SSD, 
low but significant raggedness index and negative Fu's 
Е; and Tajima's D values indicate these populations 
as a whole have possibly undergone recent expansion. 
The Bayesian skyline plot of E. argus suggested that 
population size remained relatively stable or grew slowly 
through time during both cold (positive 6"O values) and 
warm periods (negative 6"O values) (Figure 5). Increase 
was followed by a decline starting at about 0.0496 Ma, 
and then decline was followed by an increase starting at 
0.004 39 Ma. 

For all populations of Е. brenchleyi, negative Ға” 
Е; values indicate populations as a whole have possibly 
undergone recent expansion. The mismatch distribution 
showed multimodal, and the P-value for SSD and 
the raggedness index with low-value both were non- 
significant. To construct Bayesian skyline plots (BSP) 
of E. brenchleyi, the rate of evolution of ND4 was set 
to 0.064 5 substitutions per site per Ma (from the ucld. 
mean of BSP of E. argus). The Bayesian skyline plot 
of E. brenchleyi implied that population size remained 
relatively stable or descend slowly through time during 
both cold (positive 8'*O values) and warm periods 
(negative 8'"О values) (Figure 6). Decline was followed 
by an increase starting at about 0.030 6 Ma. The TMRCA 
of E. argus and Е. brenchleyi were estimated about 0.493 
Ma (with the 95% credibility interval ranging from 0.438 
to 0.548 Ma) and 1.514 Ma (with the 95% credibility 
interval ranging from 0.922 to 2.248 Ma), respectively. 

Results of AMOVA indicated the presence of 
geographic genetic structure within both species, 
concordant with results of each BSP above. For E. 
argus, a two-level AMOVA to test for structure among 
populations detected high genetic structure (Фо = 0.594, 
P « 0.001). A small portion of the total genetic structure 
exists within populations, and the majority was among 
populations (Ф; = 0.594, P < 0.001). For Е. brenchleyi, a 
two-level AMOVA to test for structure among populations 
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Figure 5 Pleistocene ice distribution (bottom) and Bayesian skyline 
plot (BSP) (top) representing the historical demographic trend of 
mitochondrial lineages of Е. argus. The solid line represents the 
mean of the effective population size, while dotted lines represent 
the 95% HPD limits. The 8:0 curve is used as proxy for global ice 
distribution (data from the SPECMAP project: http://www.ngdc. 
noaa.gov/mgg/geology/specmap.html). Negative values represent 
less ice; positive values represent more extensive ice coverage. 


detected relatively higher genetic structure (Фо = 0.908, 
P « 0.001). A small portion of the total genetic structure 
exists within populations, and the majority was among 
populations (®,; = 0.905, P < 0.001). 

For E. argus, pairwise area estimates of Ф; ranged 
from 0.149 to 0.868. АП of the 36 pairwise Фо, values 
were statistically significant (Table 3). Results of the 
Mantel test between Фо and geographical distance 
showed a weak relationship between population structure 
and distance (R^ — 0.18, P — 0.096). For E. brenchleyi, 
pairwise area estimates of ®,, ranged from 0.525 to 
0.972. АП of the 10 pairwise Фо values were significant 
(Table 4). The results of the Mantel test between Фо and 
geographical distance showed a significant relationship 
between population structure and distance (R^ — 0.43, 
P = 0.007). Thus, only Е. brenchleyi shows conclusive 
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Figure 6 Pleistocene ice distribution (bottom) and Bayesian 
skyline plot (top) representing the historical demographic trend of 
mitochondrial lineages of Е. brenchleyi. The solid line represents the 
mean of the effective population size, while dotted lines represent 
the 95% HPD limits. The oxygen isotopic curve (8"O) is plotted as 
proxy for global ice distribution (data from the SPECMAP project: 
http://www.ngdc.noaa.gov/mgg/geology/specmap.html). Negative 
values represent less extensive ice coverage; positive values 
represent more extensive ice coverage. 


evidence of isolation by distance. 

Table 5 shows results of SAMOVA. For E. argus group 
structure was maximized at k = 7 (Фет = 0.504). Тһе 
maximum among-group structure was: JY and CA vs. 
HRB vs. XH and YQ vs. GH vs. ETK vs. HD vs. CZ. For 
E. brenchleyi, there were distinct groups of genetically 
defined sampling areas, with c, ranging from 0.489 to 
0.672 when the program was instructed to identify k — 2 
through & = 4. Partitions of sampling areas were identified 
that suggested SYR vs. NYR groups (partitions: YO, HD 
and JY vs. TS and XZ; Ф.т = 0.708) at k = 2, and this 
grouping pattern was concordant with the division of the 
Yellow River and Bayesian phylogenetic tree (Figure 3). 
Group structure was maximized at k = 3 (Ф.т = 0.808). 
However, partitions of sampling areas suggest YO vs. HD 
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Table 3 Eremias argus pairwise estimates of genetic structure (Ф;т), + represents significance at P < 0.05. 


HRB XH ETK YQ HD GH JY CA (22, 
HRB 0 
XH 0.702 0 
ETK 0.8104 0.8354 0 
YQ 0.673 0.3344 0.7867 0 
HD 0.6144 0.377 0.7324 0.4074 0 
GH 0.726 0.798 0.868 0.598 0.628 0 
JY 0.6664 0.5824 0.792 0.5394 0.2854 0.767-- 0 
СА 0.521-- 0.382-- 0.6364 0.4264 0.2294 0.5454 0.1494 0 
CZ 0.5564 0.489+ 0.706+ 0.470+ 0.329+ 0.604+ 0.369+ 0.310+ 0 


Table 4 Eremias brenchleyi pairwise estimates of genetic structure 
(Фа), + represents significance at P < 0.05. 


YQ HD TS JY XZ 
YQ 0 

HD 0.525+ 0 

TS 0.940- 0.879+ 0 

ТҮ 0.723+ 0.531+ 0.970- 0 

х2 0.956- | 0.922- 0.960+ 0.9724 0 


and JY vs. TS vs. XZ [respectively corresponding to clade 
YQ, (clades HD1, HD2, HD3, HD4 and JY), clade TS 
and clade XZ in the Bayesian phylogenetic tree (Figure 
3)] at k = 4, with only a minimal reduction in group 
structure (Ф.т = 0.769). 


4. Discussion 


4.1 Phylogeographical patterns Bayesian analysis 
indicated that E. argus was composed of eight matrilines 
and had no geographical specificity. Unlike Е. argus, 
there were two reciprocally well-supported lineages 
(NYR and SYR) within E. brenchleyi, and their 
geographical separation occurred at the Yellow River 
(Figure 1, Figure 3). Results of MJN analysis revealed 
geographic genetic structuring in the two species, 
coinciding with results of the Bayesian analysis. There are 
no evident signs of isolation in Е. argus, but the Yellow 
River and Taihang Mts may have acted as important 
barriers to gene flow in E. brenchleyi. SAMOVA analysis 
of E. brenchleyi indicated that overall haplotypes were 
clustered into four distinct areas, corresponding to the 
four phylogeographical units (Avise, 2000) and and 
concordant with the Bayesian and MJN analyses. 

The Yellow River formed during the initial stage of 
the early Pleistocene (about 1.6 Ma) (Li et al., 1996; 
Zhang et al., 2003), and the divergence time between 
the two clades in E. brenchleyi was estimated at 3.208 + 
0.743 Ма (with the 95% credibility interval ranging from 
1.648 to 4.309 Ma). Accordingly, we conclude that Е. 


brenchleyi may have been separated by the Yellow River 
into two lineages (NYR and SYR). This conclusion is 
consistent with that drawn previously using cyt 5 gene 
(Zhao et al., 2011). Within clade NYR, the divergence 
time of the (HD1+HD2+HD3+HD4+JY) group from 
the YQ group was dated to 0.939 + 0.418 Ma (with the 
95% credibility interval ranging from 0.283-1.884 Ма). 
Therefore, the NYR lineage may have been separated 
by the Taihang Mts into two groups (Figure 1, Figure 
3) because Taihang Mts was uplifted to 1100-1500 m 
during the Pleistocene-Holocene (Wu et al., 1999). The 
SYR lineage may have been separated into two groups 
(TS and CZ), although we did not detect the mountain or 
river with which the geographical division of the clades 
TS and CZ coincided. The mantel test revealed that Е. 
brenchleyi had more strongly restricted gene flow among 


populations than did Е. argus. We suppose that genetic 
divergence results primarily from the lack of gene flow 
between the two regions (TS and CZ), probably because 
they are far enough from each other to result in isolation 
by distance. Dispersal ability is greater in Е. argus than 
in E. brenchleyi (Zhao et al., 2011), so the Yellow River 
may be less likely to serve as a barrier for Е. argus to 
disperse between SYR and NYR during Yellow River 
zeroflow, and the Taihang Mts are also less likely to act 
as a boundary between the lineages of E. argus. This 
conclusion is consistent with that drawn for Е. argus 
based оп суі b gene data (Zhao et al., 2011). 

АМОУА analysis indicated that among-population 
structure was high in Е. argus and even higher in Ё. 
brenchleyi, consistent with many studies showing that 
geographic genetic structuring is strong for mtDNA 
variation in lizards (Clark et al., 1999; Stenson 
and Thorpe, 2003; Jin et al., 2008; Urquhart et al., 
2009). The range of genetic structure is narrower in 
E. argus (from 0.149 to 0.868) than in E. brenchleyi 
(from 0.525 to 0.972). Eremias brenchleyi had more 
significant geographical structuring of haplotype diversity 
than Е. argus. 
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Table 5 Result of SAMOVA using mitochondrial ND4 sequences. The best number of groups has the largest Por value 


Amon Among nen Percentage of variation (96) 
Speci G G iti rou 5 populations о uum Among ithi 
pecies тоирѕ тоир composition group within groups pop Among populations 2. 
Ds Фо EIOUPS within groups popuraton 
E. argus ы OI ine 2. GĦ, - sol 0.519*** — 0.76*** 50.13 25.89 23.98 
XH, YQ, HD, GH, JY, СА, and * M m 
Three groups CZ vs. HRB vs. ЕТК 0.436 0.445 0.687 43.56 25.1 31.34 
HD, JY, CA, and CZ vs. XH, YQ, EA к Ші 
Four groups and GH vs. ETK vs. HRB 0.432 0.356 0.634 43.16 20.22 36.62 
қ HD, JY, and СА vs. НЕВ vs. ЕТК ій xd — 
Five groups vs. CZ vs. XH, YO, and GH 0.449 0.311 0.62 44.89 17.14 37.97 
M ETK vs. XH and YQ vs. HD, JY, AK — ok 
Six groups and CA vs. HRB vs. CZ vs. GH 0.484 0.253 0.615 48.38 13.08 38.55 
JY and CA vs. HRB vs. XH and xi ak UR 
Seven groups YQ vs. GH vs. ETK vs. HD vs. CZ 0.504 0.203 0.604 50.39 10.05 39.55 
, - XH vs. ETK vs. YQ vs. HD vs. GH " xx — 
Eight groups vs. CZ vs. JY and CA vs. HRB 0.499 0.202 0.6 49.89 10.14 39.97 
E. brenchleyi Two groups YQ, HD, and JY vs. TS and XZ 0.708 0.766*** 0.932*** 70,84 22.35 6.81 
Three groups TS vs. YQ, HD, and JY vs. XZ 0.808 0.65*** 0.933*** 80.77 12.49 6.74 
Four groups YQ vs. HD and JY vs. TS vs. XZ 0.769 (0:615 *** 0.911*** 76,87 14.22 8.91 


* indicates significance at P « 0.05; ** indicates significance at P « 0.01; *** indicates significance at P « 0.001. 


4.2 Genetic diversity and multiple refugia Our data 
show that neither in Е. argus nor E. brenchleyi are 
diversity indices correlated with latitude and that genetic 
diversity is almost homogeneous across the range of 
each species. These findings are consistent with cyt b 
data reported previously for the two species (Zhao et al., 
2011) but not with phylogeographic studies of widespread 
terrestrial species in North America and Europe, where 
they expanded from southern refugia at the end of 
the Pleistocene, thus exhibiting more geographically 
structured genetic diversity in the south (Hewitt, 1996, 
1999, 2000). It may be the case that multiple refugia 
were maintained across the range of each species during 
the LGM, with neither species expanding from a single 
refuge. 

Bayesian analysis shows that the haplotype clade of E. 
argus is composed of six major clades (EAIIal, EAIIa2, 
EAIIb, EAIIc, EAIId and EAIIe) with strong nodal 
support (all PP > 0.7). The TMRCA of these six clades 
is estimated to be 0.158 + 0.063, 0.178 + 0.064, 0.205 + 
0.087, 0.220 + 0.082, 0.144 + 0.082 and 0.184 + 0.082 
Ma, respectively. The haplotype clade of E. brenchleyi 1s 
composed of five major clades (JY, YO, HD4, TS and XZ) 
with strong nodal support (all PP > 0.8). The TMRCA of 
these five clades is estimated to be about 0.735 + 0.394, 
0.757 + 0.369, 0.575 + 0.355, 0.404 + 0.324 and 0.262 + 
0.264 Ma, respectively. Thus, all major clades for these 
two lizards formed before the LGM (about 0.06—0.1 Ma) 
(Tian et al., 2009). Considering that ancestral haplotypes 


are more likely to be found at refugia (Slatkin , 1996; 
Parsons et al., 1997; Li et al., 2009), and the TMRCA 
of these clades, it is possible that six major clades of Е. 
argus and five major clades of Е. brenchleyi diagnose 
six separate glacial refugia in E. argus and five separate 
glacial refugia in E. brenchleyi during the LGM. 
Furthermore, we propose an alternative hypothesis 
that each species expanded from a single refuge after the 
LGM. In this case, related haplotypes could be expected 
to form a star-like network circling the ancestral haplotype 
(Avise, 2000), and among-population structure could be 
expected to be low (Hewitt, 2000). In contrast to these 
expectations, we detected a reticulated rather than star- 
like network in both species (Figure 4), and found high 
among-population structure in Е. argus (Фо, = 0.594, P < 
0.001) and even higher among-population structure in Е. 
brenchleyi (s, = 0.908, Р < 0.001) that probably resulted 
from genetic drift favouring/fixing different alleles among 
multiple independent refugia during the LGM (Tian 
et al., 2009). Multiple isolated refugia promote genetic 
differentiation between current populations (Avise, 2000). 
Our findings support the existence of multiple refugia in 
both E. argus and Е. brenchleyi, and are consistent with 
the results reported in similar phylogeographic studies 
(Hewitt, 1996, 2000; Zamudio and Savage, 2003; Li et al., 
2009; Tian et al., 2009). Multiple glacial refugia may 
have been available to East Asian species throughout 
their entire ranges (Li et al., 2009), probably because East 
Asia potentially hosts microclimatic zones capable of 
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supporting a variety of habitats in relative stability (Qian 
and Ricklefs, 2000). 


4.3 Demographic history Populations of Е. argus and 
E. brenchleyi showed genetic imprints of demographic 
expansion as inferred by a negative Fs test, which was 
also supported by the multimodal mismatch distributions 
analysis, Rogers test of sudden expansion model and 
AMOVA (Rogers et al., 1996; Castoe et al., 2007; Nuñez 
et al., 2011). According to the inference of population 
expansion, Е. argus has undergone historical slow growth 
or is stable in size and Е. brenchleyi has undergone 
historical slow decline in size. 

АП populations of Е. argus remained relatively stable 
or grew slowly during cold (positive 5'°O values) and 
warm (negative 5'°O values) periods, probably because 
the species survived in multiple refugia during cold 
periods (stadials). Occasional migration among expanded 
populations might have resumed during warmer periods 
(interstadials), and this intermittent gene flow between 
refugia might have made the species population size 
insensitive to population reduction resulting from cool 
climates (Li et al., 2009). 

All populations of E. brenchleyi remained relatively 
stable or descended slowly through time during cold and 
warm periods, probably because there had been lower 
intermittent gene flow among populations in Е. brenchleyi 
than in Е. argus. Moreover, the expansion events 
beginning at 0.004 39 Ma for E. argus and 0.030 6 Ma 
for E. brenchleyi suggest that climate changes during the 
LGM may have differentially affected the demographic 
histories of these two species. 
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Appendix А Eremias argus haplotype frequency by population. Appendix B Eremias brenchleyi haplotype frequency by 


Haplotype labels follow Figure 2 and Figure 4 A. population. Haplotype labels follow Figure 3 and Figure 4 B. 
Haplotype Haplotype 

fauinber HRB XH ETK YQ HD GH JY CA CZ Humber YQ HD TS JY XZ 
ЕА1 1 ЕВІ 3 

ЕА2 2 ЕВ2 1 

EA3 4 EB3 1 

ЕА4 1 ЕВ4 2 

ЕА5 1 ЕВ5 1 

ЕА6 1 ЕВ6 1 

ЕА7 1 EB7 1 

ЕА8 1 EB8 1 

ЕА9 1 ЕВ9 3 

ЕА10 1 15 ЕВ10 1 

ЕА1 1 EBII 1 

EA12 1 ЕВ12 1 

EA13 3 EB13 1 

EA14 1 ЕВ14 1 

ЕА15 2 ЕВ15 1 

EA16 2 ЕВ16 2 

ЕА17 2 EBI7 1 

ЕА18 3 ЕВ18 1 

ЕА19 1 ЕВ19 1 

ЕА20 3 ЕВ20 1 

ЕА21 1 ЕВ21 2 

ЕА22 1 ЕВ22 1 

ЕА23 1 4 ЕВ23 2 

EA24 1 EB24 1 

EA25 1 ЕВ25 1 

ЕА26 1 ЕВ26 1 

ЕА27 1 ЕВ27 2 
ЕА28 5 ЕВ28 3 
ЕА29 5 ЕВ29 3 
EA30 2 EB30 1 
EA31 3 EB31 1 
EA32 1 EB32 1 
EA33 2 EB33 1 
EA34 1 6 

ЕА35 1 6 

EA36 4 

EA37 1 

ЕА38 1 

EA39 2 

ЕА40 

ЕА41 15 2 

ЕА42 1 

ЕА43 1 


ЕА44 f 


Appendix C GenBank accession numbers for the two lizard 
species in this study. GAN: GenBank accession number. 


Haplotype Haplo 

Mid GAN xh el T GAN 
EAI JN561171 EBI JN561215 
EA2 JN561172 EB2 JN561216 
EA3 JN561173 EB3 JN561217 
EA4 JN561174 EB4 JN561218 
EAS JN561175 ЕВ5 JN561219 
EA6 JN561176 EB6 JN561220 
ЕА7 JN561177 EB7 JN561221 
ЕА8 JN561178 ЕВ8 JN561222 
EA9 JN561179 EB9 JN561223 
ЕА10 JN561180 EBIO JN561224 
EAII JN561181 EBII JN561225 
ЕА12 7561182 EBI2 JN561226 
EA13 JN561183 EBI13 JN561227 
ЕА14 JN561184 EB14 JN561228 
ЕА15 JN561185 ЕВ15 JN561229 
ЕА16 JN561186 EBI6 JN561230 
ЕА17 JN561187 EBI7 JN561231 
ЕА18 JN561188 EB18 JN561232 
EAI9 JN561189 EBI9 JN561233 
EA20 JN561190 EB20 JN561234 
EA21 JN561191 EB21 JN561235 
EA22 JN561192 EB22 JN561236 
EA23 JN561193 EB23 JN561237 
EA24 JN561194 EB24 JN561238 
EA25 JN561195 EB25 JN561239 
EA26 JN561196 EB26 JN561240 
EA27 JN561197 EB27 JN561241 
EA28 JN561198 EB28 JN561242 
EA29 JN561199 EB29 JN561243 
EA30 JN561200 EB30 JN561244 
EA3I JN561201 EB31 JN561245 
EA32 JN561202 EB32 JN561246 
EA33 JN561203 EB33 JN561247 
EA34 JN561204 

EA35 JN561205 

EA36 JN561206 

EA37 JN561207 

EA38 JN561208 

EA39 JN561209 

EA40 JN561210 

ЕА41 JN561211 

EA42 JN561212 

EA43 JN561213 


EA44 JN561214 


